Mathematical operations and equation solving with reconfigurable metadevices

Performing analog computations with metastructures is an emerging wave-based paradigm for solving mathematical problems. For such devices, one major challenge is their reconfigurability, especially without the need for a priori mathematical computations or computationally-intensive optimization. Their equation-solving capabilities are applied only to matrices with special spectral (eigenvalue) distribution. Here we report the theory and design of wave-based metastructures using tunable elements capable of solving integral/differential equations in a fully-reconfigurable fashion. We consider two architectures: the Miller architecture, which requires the singular-value decomposition, and an alternative intuitive direct-complex-matrix (DCM) architecture introduced here, which does not require a priori mathematical decomposition. As examples, we demonstrate, using system-level simulation tools, the solutions of integral and differential equations. We then expand the matrix inverting capabilities of both architectures toward evaluating the generalized Moore–Penrose matrix inversion. Therefore, we provide evidence that metadevices can implement generalized matrix inversions and act as the basis for the gradient descent method for solutions to a wide variety of problems. Finally, a general upper bound of the solution convergence time reveals the rich potential that such metadevices can offer for stationary iterative schemes.

The basic module that is used throughout the analysis is the MZI constructed as shown in Fig. (S1). This module consists of two 3 dB directional couplers and two phase shifters. The MZIs can be assembled and configured in a variety of ways to perform matrix multiplication, i.e. the Miller and DCM architectures. Although in our simulations we used the 1 GHz frequency, this choice is arbitrary and obviously, depending on the MZI platforms, other frequencies can be used.
As seen in Fig. (S1) the module was simulated as a combination of two 3 dB couplers (QHYB) and two phase shifters (PHASE), all centered to operate at 1 GHz. The QHYB coupling module is a typical symmetrical directional coupler as the one introduced in 2 .
In the case of the multiplier module we used an MZI module followed by an amplifier, as can be seen in Fig. (S2). This module is the main building block for the DCM architecture.

Feedback coupler
Equation solving is enabled by routing the output of the matrix multiplication stage back toward the input (i.e. feedback). However, a system input must be provided and the output must be probed, all while perturbing the feedback as little as Figure S1. Basic MZI module used for the system simulation. The module was simulated as a combination of two 3 dB couplers (QHYB) and two phase shifters (PHASE), all centered to operate at 1 GHz. Notice that the QHYB coupling module is a typical symmetrical directional coupler as the one introduced in 2 .
possible. To achieve this, we use pairs of the four-port coupler presented in Fig. (S3). The coupler implements the following where α and β are the transmission and coupling coefficients, respectively. In all examples |β | 2 = −20 dB, i.e., 1% of the power is coupled to the output port while the rest (|α| 2 = 1 − |β | 2 = 99%) is directed to the system, while for the inversion of singular and rectangular cases the coupling coefficient was |β | 2 = −30 dB (0.1%). Notice that in both cases we utilized two couplers, allowing extra control over the extracted results and isolation of the output signal from the input signal, similar to 3 .

Miller architecture
The basic modules can be combined to form the Miller architecture to perform matrix multiplication. With the addition of the feedback loop and coupling elements, linear equation solving can be achieved. This is depicted in Fig. (S4) for a 7 × 7 discretized system. Note that this architecture comprises of 28+7+28=63 MZI modules. In this particular example we utilized two different sets of feedback couplers, used for system inputs and outputs, respectively. Both sets had identical coupling/transmission coefficients. Figure S2. Basic MZI module followed by an amplification unit representing the multiplier module.

The DCM architecture
An alternative to the Miller architecture for matrix multiplication is the DCM architecture. Similar with the above, Fig. (S4) depicts the 5 × 5 equation solving system used to extract the presented results in the main text. For this system, a total of 5 2 = 25 MZIs are used for the middle stage. Both the ingress and egress stages (Fig. (S5)) comprise a total of 5 × 4 MZIs with fixed phase shifts ratios. These devices are passive and are necessarily lossy. Both ingress and egress stages introduce a signal division of the order of α in = α en = 1/ √ 5. The power losses associated with the ingress and egress stages are −10 log(5 2 ), and can be compensated with an amplification stage (that is part of the main MZI module) that introduces a total of 10 log(5 2 ) = 13.97 amplification to the system. Additionally, the feedback coupler introduces a α fc loss to the system. Therefore, the DCM system can be described by the following equations and in this case the resulted inversion will yield to the following solution In the ideal case the transmission coefficients α in α eg α fc can all be compensated by the amplification stage with gain g, hence the overall solution will be x ≈ A −1 e in , assuming that the overall required amplification level is g = 10 log(5 2 /|α fc | 2 ). Figure S3. Module that is used for the feedback couplers, centered to operate at 1 GHz.

The ingress and egress stages
For the particular example of an 11 × 11 network, both the ingress and egress stages were simulated as a cascade of 10 MZI with fixed splitting ratios, i.e., the 1st MZI had a 1:10 balance, the 2nd 1:9, ..., and the 10th had 1:1 balance, designed so that if the cascade is illuminated with 1W of power, one would find 1/11 W on the outputs of each splitter. Both of these stages can be seen in Fig. (S6)

Results retrieval
From a system point of view, the solution (regardless of the architectures) can be retrieved by simply probing the input channels, just before the main entrance operator. In reality such measurement will disturb the actual result. For these reasons, the results are extracted using a second set of feedback couplers, as seen in Fig. (S4) and Fig. (S5), top-left corner. Assuming that the second set of couplers have the same coupling/transmission ratios as the first, then the extraction of the results is a straightforward task, similar to the one presented in 3 . Using a separate set of couplers for the input and output prevents the uncoupled input energy from overwhelming the output signal. Figure S4. A schematic of the Miller architecture, centered to operate at 1 GHz. Note that in this scheme two feedback couplers were used, one for the input and one for the output of the system. Figure S5. A schematic of the 5 × 5 DCM architecture all centered to operate at 1 GHz. Figure S6. A detail of the the 5 × 5 DCM architecture depicting the the ingress, middle, and egress stages.

Statistical comparison of Miller and DCM architectures
Here we include an example where we performed a matrix inversion of a rectangular complex-valued 3x5 matrix with a series of tolerance/yield analyses for both architectures. In particular, we examine what is happening when we assume that the MZIs (phase-shifting elements) exhibit a normal distribution variation applied to the phase shifter and to the specific phase value.
Note that for DCM architecture, both ingress and egress stages could be implemented using variable MZI components, but since the coupling ratios are fixed, we are assuming that these were implemented using fixed perfect splitters, hence we exclude them from our statistical analysis. Additionally, we did not include the effects of the amplifier noise, as this would confound the comparison between the two architectures.
First, let's examine the new complex-valued rectangular matrix. The randomly chosen matrix reads For the matrix inversion, we use the same system-level simulation process used in the text. The numerical pseudoinverse of this matrix reads In particular, we have assumed that the kernel is with α λ = 0.25. Due to the rectangular nature of the matrix A, the smallest eigenvalues of A * A are zero. Therefore, when performing the SVD of I − α λ A * A, as required in the Miller architecture, we find that the corresponding two larger singular values of the kernel are s 1 = s 2 ≈ 1. Having eigenvalues that are very close to unity makes the convergence of the system very poor, yielding noise sensitivity that may lead to inaccurate results. For these reasons in the Miller architecture, since these two singular values can affect its convergence, we reduce these values to be s 1 = s 2 = 0.5. In other words, since we know the 7/11 spectral distribution of the system (due to the SVD) we tune the extreme eigenvalues of the kernel to be within the unit circle.
As we will see, not only does this choice improve the accuracy of the results for the Miller architecture, but also it endows the Miller architecture with robust characteristics.
Next, we endeavor to reduce the sensitivity of the DCM. Obviously, we cannot do the same for DCM as was done for the Miller architecture since we do not have access to any kind of spectral information of the matrix -no a-priori mathematical operation is done to the kernel. Instead, we introduce a small perturbation to the kernel, i.e, where δ is a small value -these kinds of perturbations are typically encountered in Tikhonov-type regularization schemes (see Methods in text). In this sense, we reduce the sensitivity of DCM by introducing small reduction. Since for the case of rectangular/singular matrix A, one or more of the eigenvalues is zero, this reduction will pull the eigenvalues of I − α λ A * A slightly away from the rim of the unit circle. The results indicate that the Miller architecture has a better response to phase variations when compared to DCM. This result demonstrates that Miller's architecture is more robust and tolerant to noise. These features are mainly due to the SVD decomposition. DCM is slightly more sensitive to phase variations -in DCM, we need to use a gain stage that also amplifies the various effects of the overall system.